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1 Introduction 

Theoretical study of polarons in the strongly correlated system is like an at- 
tempt to view contents of a Pandora box embedded into another, even more 
sinister and obscure, container of riddles, enigmas and mysteries. This des- 
perate situation occurs because solution is not known even for the simplest 
polaron problem, i.e. when a perfectly stable quasiparticle (QP) with mo- 
mentum as a single quantum number interacts with a well defined bath of 
bosonic elementary excitations. To the contrary, the definition of the strongly 
correlated system implies that QPs might be highly unstable and the very 
notion of QPs, both in electronic and bosonic subsystems, is under question. 
Thus, one faces the problem of an interplay between ill defined objects and 
it is crucial to solve the problem without approximations. Further difficulty, 
pertinent to realistic systems, is an interplay of the momentum and other 
quantum numbers characterizing internal states of a QP. 

The problem of polaron originally emerged as that of an electron coupled 
to phonons (see [1, 2]). In the initial formulation a structureless QP is char- 
acterized by the only quantum number, momentum, which changes due to 
interaction of the QP with phonons [3, 4]. Later, depending on what can be 
called "particle" and "environment", and how they interact with each other, 
the polaron concept was related to extreme diversity of physical phenomena. 
There are many other objects which, having nothing to do with phonons, are 
isomorphic to simple polaron [5], as, e.g. an exciton-polaron in the intraband 
scattering approximation [6, 7, 8, 9]. Another example is the problem of a hole 
in the antiferromagnet which is closely related to polaron since hole movement 
is accompanied by the spin flips which, in the spin wave approximation, are 
equivalent to creation and annihilation of magnons [10, 11]. 
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The concept of polaron was further generalized to include internal degrees 
of freedom which, interacting with environment, change their quantum num- 
bers. Example of a complex QP is the Jahn- Teller polaron, where electron- 
phonon interaction (EPI) changes quantum numbers of degenerate electronic 
states [12, 13, 14]. This generalization is important due to it's relevance to 
the colossal magnetoresistance phenomena in the manganese oxides [15, 16]. 
Another example is the pseudo Jahn- Teller polaron, where EPI is inelastic 
and leads to transitions between close in energy electronic levels of a QP 
[17, 18, 19]. Further generalization is a system of several QPs which interact 
both with each other and environment. For example, effective interaction of 
two electrons through exchange by phonons can overcome the Coulomb repul- 
sion and form a bound state, bipolaron [20, 21, 22, 23, 24]. On the other hand, 
coupling of attracting hole and electron to the lattice vibrations [25, 26, 27] 
can create a lot of qualitatively different objects: localized exciton, weakly 
bound pair of localized hole and localized electron, etc. [28, 7]. Scattering by 
impurities introduces additional complexity to the polaron problem because 
interference of impurity potential with lattice distortion, which accompanies 
the polaron movement, can contribute cither constructively or destructively 
to the localization of a QP on impurity [29, 30, 7]. 

In addition, a bare QP and bosonic bath can not be considered as well 
defined in the correlated systems. Angle Resolved Photoemission Spectra 
(ARPES), revealing the Lehmann Function (LF) of quasiparticle, demonstrate 
broad peaks in many correlated systems: cooper oxide high-temperature su- 
perconductors [31, 32, 33], colossal magnetoresistivc manganites [34, 35, 36], 
quasi-one-dimensional Peierls conductors [37, 38] , and Verwey magnetites [39] . 
Besides, phonons are also broadened in many correlated systems, e.g. in high- 
temperature semiconductors [40] and mixed- valent materials [41, 42]. One of 
possible reasons for these broadenings is the interaction of the QPs with the 
lattice degrees of freedom. However, in many realistic cases other subsystems, 
not explicitly included into the polaron Hamiltonian, are responsible for the 
decay of QP and phonons, e.g., another electronic bands, phonon anharmonic- 
ity, interaction with nuclear spins, etc. Then, if this auxiliary broadening is 
known in some approximation, one can formulate an ambitious goal to study 
spectral response when "bare" quasiparticle with known damping interacts 
with "broadened" bosonic excitations. 

No one of traditional numerical methods, to say nothing of analytical ones, 
can give approximation free results for measurable quantities of polaron, such 
as optical conductivity or angle resolved photoemission spectra, for in macro- 
scopic system of arbitrary dimension. Besides, we are not aware of any numer- 
ical method which can incorporate in an approximation free way the informa- 
tion on the damping of QP and bosonic bath. Below we describe basics of re- 
cently developed Diagrammatic Monte Carlo (DMC) method for numerically 
exact computation of Green functions and correlation functions in imaginary 
time for few polarons in a macroscopic system [43, 44, 45, 46, 47, 48, 49, 50, 51]. 
Analytic continuation of imaginary time functions to real frequencies is per- 
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formed by a novel approximation free approach of stochastic optimization 
(SO) [45, 50, 51], circumventing difficulties of popular Maximal Enthropy 
method. Finally we focus on results of application of the DMC-SO machinery 
to various problems [52, 53, 54, 55, 56, 57] 

The basic models, related to the polaronic objects in correlated systems, 
which can be solved by DMC-SO methods, are stated in the next Sect. It 
is followed in Sect. 1.2 by description of stumbling blocks encountered by 
analytic methods. Sect. 2 concerns the basics of DMC-SO methods. However, 
those who are not interested in the details of the methods can briefly look 
through the definitions in the introduction for Sect. 2 and turn to Sect. 3 
where LF and optical conductivity of Frohlich polaron are discussed (see also 
[58]). Results of studies of the self-trapping phenomenon are presented in Sect. 
4 and application of DMC-SO methods to the exciton problem can be found 
in Sect. 5. The chapter is completed by Sect. 6 devoted to studies of ARPES 
of high temperature superconductors. 



1.1 Formulation of a General Model with Interacting Polarons 

In general terms, the simplest problem of a complex polaronic object, where 
center-of-mass motion does not separate from the rest of degrees of freedom, 
is introduced as system of two QPs 

^ par = E e «( k ) a k«k + 5> h (k)fc k 4 (i) 

k k 

(a k and /ik are annihilation operators, and e a (k) and £/j(k) are dispersions of 
QPs), which interact with each other 

i?a-h - -A^ 1 ]T W(p,k,k')4 +k ^_ k V-k<a P +k'. (2) 
P kk' 

(N is the number of lattice sites) through the instantaneous Coulomb potential 
and the scattering by bosons 

Q 

K=l k,q 

7oa, K (k,q)aj c _ q a k + 7^ ire (k,q)/ij c _ q /i k + 7 o/liK (k,q)/ij c _ q a k + h.c. (3) 

(l[aa.ah,hh],K are interaction constants) where quanta of Q different branches 
of bosonic excitations are created or annihilated, which are described by 

Q 

H bos = ^Kb^ K b^ K . (4) 

re=l q 



In general, each QP can be a composite one with internal degree of freedom 
represented by T different states 
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C=^f fl (k)4 kttl , k , (5) 

k i=l 

which quantum numbers can be also changed due to nondiagonal part of 
particle-boson interaction 

Q T 

iWbos = ^EE 7y>(k, q)(& q ,« - fe-q,K)a!,k-q a J:k + h.c. (6) 

K k.q i,j = l 

Complicated model (l)-(6) is still too far from the cases encountered in 
strongly correlated systems. Due to coupling of QPs (1) and (5) and bosonic 
fields (4) to additional degrees of freedom, these excitations are not well de- 
fined from the onset. Namely, the dispersion relation of the QP spectrum e(k) 
in realistic system is ill-defined. One can speak of a Lehmann Function (LF) 
[59, 60, 61] of a QP 

Lt(w) = - E ^)) IH4|vac)| 2 (7) 

V 

, which is normalized to unity / +O ° dioL-^{uj) = 1 and can be interpreted as 
a probability that a QP has momentum k and energy w. (Here {\v}} is a 
complete set of eigenstates of Hamiltonian H in a sector of given momentum 
k: H |^(k)) = E v (k) |i/(k)).) Only for noninteracting system the LF reduces 
to delta function L^ ONINT (a>) = 5(w — e(k)) and, thus, sets up dispersion 
relation w = e(k). 

Specific cases of model (l)-(6) describe enormous variety of physical prob- 
lems. Hamiltonians (1) and (2), in case of attractive potential U(p, k, k') > 0, 
describe an exciton with static screening [62, 63]. Besides, expressions (l)-(4) 
describe bipolaron for repulsive interaction [20, 21, 22, 23, 24] U(p, k, k') < 
and exciton-polaron otherwise [25, 26, 27]. The simplest model for exciton- 
phonon interaction, when only two (T = 2) lowest states of relative electron- 
hole motion are relevant (e.g. in one-dimensional charge-transfer exciton 
[64, 65, 66]), is defined by Hamiltonians (4)-(6)). The same relations (4)-(6) 
describe the problems of Jahn- Teller [all ej in Hamiltonian (5) are the same] 
and pscudo Jahn- Teller polaron. The problem of a hole in an antiferromagnet 
in spin- wave approximation is expressed in terms of Hamiltonians (4)- (6) with 
Q = 1 and T = 1. When hole also interacts with phonons, one has to take into 
account one more bosonic branch and set Q — 2 in (4) and (6). Finally, the 
simplest nontrivial problem of a polaron, i.e. of a structureless QP interacting 
with one phonon branch, is described by noninteracting Hamiltonians of QP 
H pal and phonons H p h 

i? = ^ e(k)4a k + ^ q 6 q 6 q , (8) 

k q 



and interaction term 
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Hint = E V ^ 1) {b \ - 6 -q)4- q a k + h.c. . (9) 

k,q 

The simplest polaron problem, in turn, can be subdivided into continuous and 
lattice polaron models. 

1.2 Limitations of Analytic Methods in Problem of Polarons 

Analytic solution for the problem of exciton in a rigid lattice is available 
only for small radius Frenkel regime [67] and large radius Wannier regime 
[68]. However, even limits of validity for these approximations are not known. 
Random phase approximation approaches [62, 63], are capable of obtaining 
some qualitative conclusions for intermediate radius regime though its' quan- 
titative results are not reliable due to uncontrolled errors. The situation is 
similar with the problem of structureless polaron, where analytic solutions 
arc known only in the weak and strong coupling regimes. Besides, reliable 
results for these regimes are available only for ground state properties. 

Although several novel methods, capable of obtaining properties of excited 
states, were developed recently, variational coherent-states expansion [69] and 
free propagator momentum average summation [70] as a few examples, all of 
them, to provide reliable data in a specific regime, need either comparison 
with exact sum rules [71, 72] or with exact numerical results. 

Application of variational methods to study of excitations is a tricky is- 
sue since, strictly speaking, they are valid only for the ground state. As an 
example for the importance of sum rules in variational treatment, we refer 
to the problem of the optical conductivity of the Frochlich polaron. Possibil- 
ity of existence of Relaxed Excited State (RES), which is a metastable state 
where lattice deformation has adjusted to the electronic excitation rendering 
stability and narrow linewidth of the spectroscopic response, was briefly men- 
tioned by S. I. Pekar in early 50's. Then, conception of RES was rigorously 
formulated by J. T. Devreese with coworkers and has been a subject of ex- 
tensive investigations for years [5, 73, 74, 75, 76, 77, 48, 57]. Calculations of 
impedance [75] in the framework of technique [78] supported the existence of 
a narrow stable peak in the optical conductivity. However, even the authors 
of [75] were skeptical about the fact that the width of RES in the strong 
coupling regime appeared to be more narrow than the phonon frequency, i.e. 
inverse time which is, according to the Hciscnberg uncertainty principle, is 
required for the lattice readaptation. In consequent paper [77] they realized 
the importance of many-phonon processes and studied two-phonon contri- 
bution to optical conductivity. Importance of many phonon processes was 
confirmed when variational results [75] were compared with exact DMC sim- 
ulations [48]. Variational result well reproduced the position of the peak in 
exact data though failed in description of the peak width in the strong cou- 
pling regime [48]. Finally, when approach [75] was modified and several sum 
rules were accurately introduced into variational model [57], both position 



6 A. S. Mishchenko and N. Nagaosa 

and width of the peak were quantitatively reproduced. Studies [57] (see Sect. 
3.1), do not address rather philosophical question whether RES exists or not, 
though inevitably prove that, in contrast to the foregoing beliefs, there in no 
stable excited state of the Frohlich polaron in the strong coupling regime. 
Note that sometimes excited states can not be handled by analytic methods 
even for weak couplings: perturbation theory expression for LF of the Frohlich 
polaron model diverges at the phonon energy w p h [See (34) in Sect. 3.1.] and 
more elaborate treatment is necessary. 

Difficulties of semianalytic methods enhance in the intermediate coupling 
regime where results are sometimes wrong even for ground state properties. 
For example, the variatioanl approach [79], which has been considered as an 
intermediate coupling theory appeared to be valid only in the weak coupling 
limit [45]. Special interest to the methods, giving reliable information on ex- 
cited states, is triggered by the self-trapping phenomenon which occurs just 
in the intermediate coupling regime. This phenomenon is a dramatic trans- 
formation of QP properties when system parameters are slightly changed 
[3, 7, 9, 80]. In the intermediate coupling regime "trapped" QP state with 
strong lattice deformation around it and "free" state with weakly perturbed 
lattice may hybridize and resonate because of close energies at some critical 
value of electron-lattice interaction j c . It is clear that, to study self-trapping, 
one has to apply a method giving reliable information on excited states in the 
intermediate coupling regime. 

2 Diagrammatic Monte Carlo and Stochastic 
Optimization Methods 

In this section we introduce definitions of exciton-polaron properties which 
can be evaluated by DMC and SO methods. An idea of DMC approach for 
numerically exact calculation of Green functions (GFs) in imaginary times 
is presented in Sect. 2.1, and a short description of SO method, which is 
capable of making unbiased analytic continuation from imaginary times to 
real frequencies, is given in Sect. 2.2. Using combination of DMC and SO, 
one can often circumvent difficulties of analytic and traditional numerical 
methods. Therefore, a brief comparative analysis of advantages and drawbacks 
of DMC-SO machinery is given in Sect. 2.3. 

To obtain information on QPs it is necessary to calculate Matsubara GF 
in imaginary time representation and make analytic continuation to the real 
frequencies [60]. For the two-particle problem (l)-(4) the relevant quantity is 
the two-particle GF [46, 47] 

G£ p ' (r) = (vac | a k+p , (r)h^ p , (r)4_ p 4 +p | vac) . (10) 

(Here hk-p(r) = e HT h^ p e^ Hr , r > 0.) In the case of exciton-polaron, vac- 
uum state | vac) is the state with filled valence and empty conduction bands. 



Spectroscopic Properties of Polarons by Exact Monte Carlo 7 

For the bipolaron problem it is a system without particles. In the simpler case 
of a QP with two- level internal structure described by (4)- (6) the relevant 
quantity is the one-particle matrix GF [52, 47] 

Gk,ij(T) = (vac | o»,k(r)ot k | vac), i,j = 1,2. (11) 
For a structureless polaron the matrix (11) reduces to one-particle scalar GF 

Gk( T ) — ( vac I a k( r ) a k I vac ) • (12) 

Information on the response to an external weak perturbation (e.g. optical 
absorption) is contained in the current-current correlation function {Jp{r)J$) 
((3/5 are Cartesian indexes). 

Lehmann spectral representation of G k (r) [60, 61] at zero temperature 

/■OO 

G k (r) = / divLu(u>)e-" r , (13) 
Jo 

with the Lehmann function (LF) L k (w) given in (7), reveals information on 
the ground and excited states. Here {\f}} is a complete set of eigenstates of 
Hamiltonian H in a sector of given momentum k: H\u(k)) = E„(k) \v(k)). 
The LF Lk(u>) has poles (sharp peaks) on the energies of stable (metastable) 
states of particle. For example, if there is a stable state at energy E(k), the LF 
reads L k (co) = 5(uj — E(k)) + . . ., and the state with the lowest energy 
Sg. s .(k) in a sector of a given momentum k is highlighted by asymptotic 
behavior of GF 

G k (r » max [w-jj ) -> exp[-£; g . s .(k)r] , (14) 

where Z^ k ^-factor is the weight of the state. Analyzing the asymptotic behavior 
of similar n-phonon GFs [45, 52] 

Gk(n,r; qi, . . . ,q„) = _ 
(vac| 6 qn (r) • • • 6 qi (r) a p (T)al 6 qi ■ ■ ■ & qn |vac) , p = k Y,]Li 

one obtains detailed information about lowest state. For example, important 
characteristics of the lowest state wave function 

T oo 

^ g . s .(k)=]T]r J2 ^( k ;qi— t l«)4 k _ qi ..._ qn 6 qi ...6 q Jvac) (16) 

i— 1 n— qi...q n 

are partial n-phonon contribution 

T 

Z^{n)^Y. E I 0i(k;qi,...,q„) | 2 (17) 

i=l qi...q„ 

which is normalized to unity Y^!=o Z^{n) = 1, and the average number of 
phonons 
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(AO ee <# g ,.(k) \J2 b % I **B.(k)> =Y,nZ^{n) 



(18) 



in polaronic cloud. Another example is the wave function of relative electron- 
hole motion of exciton in the lowest state in the sector of given momentum 



^g.s.(k) = P^^Oak+p^k-p I vac ) • 



(19) 



The amplitudes £k p(<7-s.) of this wave function can be obtained [46] from 
asymptotic behavior of the following GF (10) 



GT P V -oo)=|& P M |2 e -%.s.(k)r _ 



(20) 



Information on the excited states is obtained by the analytic continuation 
of imaginary time GF to real frequencies which requires to solve the Fredholm 
equation Gk(i") = T [L^(lj)} (13) 



£k(w)=^ 1 [Gk(r)] 



(21) 



The equation (13) is a rather general relation between imaginary time GF/cor- 
relator and spectral properties of the system. For example, the absorption 
coefficient of light by excitonsX(w) is obtained as solution of the same equation 
[46] 



(22) 



Besides, the real part of the optical conductivity 0735 (o>) is expressed [48] in 
terms of current-current correlation function (Jp(r)Js) by relation 



apsiw) = tt^ 1 [(Jp(t)Js)} /u • 



(23) 



2.1 Diagrammatic Monte Carlo Method 

DMC Method is an algorithm which calculates GF (10)-(12) without any 
systematic errors. This algorithm is described below for the simplest case of 
structureless polaron [45], and generalizations to more complex cases can be 
found in consequent references 4 . DMC is based on the Feynman expansion of 
the Matsubara GF in imaginary time in the interaction representation 



4 Generalization of described below technique to the case of exciton (1-2) is given 
in [46] and its modification for pseudo-Jahn- Teller polaron (4-6) is developed in 
[52, 47]. Method for evaluation of current-current correlation function can be 
found in [48] and a case of a polaron interacting with two kinds of bosonic fields 
is considered in [49]. 
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Gk(r) = ( vac 



a k (T)aj c (0)exp|- f H in t(r')dr'l vac\ ;r>0 

I JO J J /con 



(24) 

Here T T is the imaginary time ordering operator, |vac) is a vacuum state with- 
out particle and phonons, H lnt is the interaction Hamiltonian in (9). Symbol of 
exponent denotes Taylor expansion which results in multiple integration over 
internal variables {t[, t 2 , . . .}. Operators are in the interaction representation 
A(t) = cxp[T(_ff par + 7J p h)]Aexp[— r(H par + H p h)]- Index "con" means that 
expansion contains only connected terms where no one integral over internal 
time variables {t[, t 2 , . . .} can be factorized. 

Vick theorem expresses matrix element of time-ordered operators as a 
sum of terms, each is a factor of matrix elements of pairs of operators, and 
expansion (24) becomes an infinite series of integrals with an ever increasing 
number of integration variables 

oo 

G kW = £ £ dx' 1 ---dx' m V^\r;{x[,...,x' m }). (25) 

m=0,2,4... £ m J 

Here index £ m stands for different Feynman diagrams (FDs) of the same 
order m. Term with m = is the GF of the noninteracting QP G^\t). 

Function T>m m \T~; {x' l7 . . . , x' m }) of any order m can be expressed as a fac- 
tor of GFs of noninteracting quasiparticle, GFs of phonons, and interaction 
vortexes V(k, q). For the simplest case of Hamiltonian system expressions 
for GFs of QP G^(t 2 — n) = cxp [— e(k)(r 2 — n)] (r 2 > n) and phonons 

Dq \r 2 — Ti) = cxp [— L0q(T2 — Ti)] (t 2 > T\) are well known. 

An important feature of the DMC method, which is distinct from the 
row of other exact numerical approaches, is the explicit possibility to include 
renormalized GFs into exact expansion without any change of the algorithm. 
For example, if a damping of QP, caused by some interactions not included 
in the Hamiltonian, is known, i.e. retarded self-energy of QP S lct (k, us) is 
available, renormalized GF 

G(°) (r) = I f° due— ImSMKu) 

W-co [u-e(k)-ReS Ict {k 7 u)} 2 + [ImS rct (k,uj} 2 

can be introduced instead of bare GF (t) . Explicit rules for evaluation of 
Dm™' do not depend on the order and topology of FD. GFs of noninteracting 
QPs (t 2 — ti) (or G^ (t 2 — ti)) with corresponding times and momenta are 

ascribed to horizontal lines and noninteracting GFs of phonon Dq\r2 — t\) 
(multiplied by the factor of corresponding vortexes V(k', q)V* (k", q)) are 
attributed to phonon propagator arch (see Fig. la). Then, Dm™' is the factor 
of all GSs. For example, expression for the weight of the second order term 
(Fig. lb) is the following 
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P 2 (r;K,r(,q}) 

?L 0) W)4°2 q 



\V(k, q)| 2 D$\t> t[)G { °\t' 1 )G ( £ ) _(t> 2 - r{)4 0) (r - r' 2 ) . (27) 




Fig. 1. (a) Typical FD contributing into expansion (25). (b) FD of the second 
order and (c) forth order. 



The DMC process is a numerical procedure which, basing on the Metropo- 
lis principle [81, 82], samples different FDs in the parameter space (r, to, £ to , {cc^}) 
and collects statistics of external variable t in a way that the result of this 
statistics converges to exact GF Gk(r). Although sampling of the internal 
parameters of one term in (25) and switch between different orders is per- 
formed within the the framework of one and the same numerical process, 
it is instructive to start with the procedure of evaluation of a specific term 



V 



(U) 



(t;{xi, 



»})■ 



Starting from a set {r; {x[, . . . ,2;^}}, an update x 



(old) 



X 



(new) 
I 



of an 



arbitrary chosen parameter is suggested. This update is accepted or rejected 
according to Metropolis principle. After many steps, altering all variables, 
statistics of external variable converges to exact dependence of the term on 



t. Suggestion for new value of parameter x 



(new) _ 



S 1 (R) is generated by 



random number R G [0, 1] with a normalized distribution function W(x{) 



in a range x 



< Xl < x 



(max) 



There are only two restrictions for this 



otherwise arbitrary function. First, new parameters x[" e "^ must not violate 
FD topology, i.e., for example, internal time t[ in Fig. lc must be in the range 



x 



(min) __ 



= 0,x 



(max) 



3]. Second, the distribution must be nonzero for the 



whole, allowed by FD topology, domain. This ergodicity property is crucial 
since it is necessary to sample the whole domain for convergence to exact 



answer. At each step, update x 



(old) 



(new) 
C l 



is accepted with probability 



1 acc - 

form 



M (if M < 1) and always otherwise. The ratio M has the following 
\r;{x' 1 ,...,x ( r W \...^'m})/W(x ( r w) ) 



M 



V 



V 



(U) 



,x 



(old) 



,x' m })/W{xf d) ) 



(28) 



For uniform distribution W = const = 



(max) 



(max) 



the probability 



of any combination of parameters is proportional to the weight function V. 
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However, for better convergence the distribution W(xY ew ) must be as close as 

possible to the actual distribution given by function T>m m \{. . . , x\ new \ ...,}). 

For sampling over FDs of all orders and topologies it is enough to introduce 
two complimentary updates. Update A transforms FD T>m m ^ (r; {x[, . . . , x' m }) 
into higher order FD ^^^(t; {xi, . . . , x' m ; q', Tg, 7-4}) with extra phonon 
arch, connecting some time points T3 and t' a by phonon propagator with mo- 
mentum q' (Fig. lc). Note that the ratio of weights 5„+2 /^m™' is n °t 
dimensionless. Dimensionless Metropolis ratio 

M _PA V ( t + r\T-A^---^' m , q',T',r"» 
PBV ( i m) (T;{x[,...,x' m })W(< l \T',T") ' 

contains normalized probability function W{q[, r', r"), which is used for gen- 
erating of new parameters 5 . Complementary update £>, removing the phonon 
propagator, uses ratio M _1 [45]. 

Note that all updates are local, i.e. do not depend on the structure of the 
whole FD. Neither rules nor CPU time, needed for update, depends on the 
FD order. DMC method does not imply any explicit truncation of FDs order 
due to finite size of computer memory. Ever for strong coupling, where typical 
number of phonon propagators N p h, contributing to result, is large, influence 
of finite size of memory is not essential. Really, according to Central Limit 
Theorem, number of phonon propagators obeys Gauss distribution centered 
at Nph with half width of the order of y/N p h [83] . Hence, if a memory for at 
least 2N p h propagators is reserved, diagram order hardly surpasses this limit. 



2.2 Stochastic Optimization Method 



The problem of inverting of integral equation (13) is an ill posed problem. 
Due to incomplete noisy information about GF Gk(r), which is known with 
statistic errors on a finite number of imaginary times in a finite range [0, r max ], 
there is infinite number of approximate solutions which reproduce GF within 
some range of deviations and the problem is to chose "the best one" . Another 
problem, which is a stumbling block for decades, is the saw tooth noise insta- 
bility. It occurs when solution is obtained by a naive method, e.g. by using 
least-squares approach for minimizing deviation measure 



D[L k (u>)] 



G k (r)-G k (r) G^{t)<1t 



(30) 



Here G^{t) is obtained from approximate LF Lk(u>) by applying of integral 
operator Gu_(t) = T Ly(io) in (13). Saw tooth instability corrupts LF in the 
ranges where actual LF is smooth. Fast fluctuations of the solution L^(w) often 



5 The factor Pa/pb depends on the probability to address add/remove processes. 
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have much larger amplitude than the value of actual LF L k (w). Standard tools 
for saw tooth noise suppression are based on the early 60-es idea of Fillips- 
Tikhonov rcgularization method [84, 85, 86, 87]. A nonlinear functional, which 
suppresses large derivatives of approximate solution L^(uj), is added to the 
linear deviation measure (30) . Most popular variant of rcgularization methods 
is the Maximal Entropy Method [61]. 

However, typical LF of a QP in a boson field consists of ^-functional peaks 
and smooth incoherent continuum with a sharp border [45, 54]. Hence, sup- 
pression of high derivatives, as a general strategy of the regularization method, 
fails. Moreover, any specific implementation of the regularization method uses 
predefined mesh in the to space, which could be absolutely unacceptable for 
the case of sharp peaks. If the actual location of a sharp peak is between 
predefined discrete points, the rest of spectral density can be distorted be- 
yond recognition. Finally, regularization Maximal Entropy approach requires 
assumption of Gauss distribution of statistic errors in Gk(r), which might be 
invalid in some cases [61]. 

Recently, a Stochastic Optimization (SO) method, which circumvents 
abovementioned difficulties, was developed [45]. The idea of the SO method 
is to generate a large enough number M of statistically independent nonreg- 
ularized solutions {L^\uj)}, s = 1, M, which deviation measures £>(») are 
smaller than some upper limit D u , depending of the statistic noise of the GF 
Gic(t). Then, using linearity of expressions (13), (30), the final solution is 
found as the average of particular solutions {L^\uj)} 

M 

L k ( W ) =M- 1 £4V). (31) 

s=l 

Particular solution (lu) is parameterized in terms of sum 

4 s) M = Ex { p t }M (32) 

t=l 

of rectangles {Pt} — {h t ,wt 7 c t } with height h t > 0, width wt > 0, and center 
d- Configuration 

C = {{P t },t = l,...,K} , (33) 

which satisfies normalization condition Y^t=i htWt — 1, defines function 
Gk(r). The procedure of generating particular solution starts from stochastic 
choice of initial configuration C* mt . Then, deviation measure is optimized by 
a randomly chosen consequence of updates until deviation is less than D u . In 
addition to updates, which do not change number of terms in the sum (32), 
there are updates which increase or decrease number K. Hence, since the 
number of elements K is not fixed, any spectral function can be reproduced 
with desired accuracy. 
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Although each particular solution (w) suffers from saw tooth noise at 
the area of smooth LF, statistical independence of each solution leads to self 
averaging of this noise in the sum (32). Note that suppression of noise happens 
without suppression of high derivatives and, hence, sharp peaks and edges 
are not smeared out in contrast to regularization approaches. Therefore, saw 
tooth noise instability is defeated without corruption of sharp peaks and edges. 
Moreover, continuous parameterization (32) does not need predefined mesh in 
w-space. Besides, since the Hilbert space of solution is sampled directly, any 
assumption about distribution of statistical errors is not necessary. 

SO method was successfully applied to restore LF of Frohlich polaron 
[45], Rashba-Pekar exciton-polaron [54], hole-polaron in t-J-model [53, 49], 
and many-particle spin system [88] . Calculation of the optical conductivity of 
polaron by SO method can be found in [48] . SO method appeared to be helpful 
in cases when GF's asymptotic limit, giving information about ground state, 
can not be reached. For example, sign fluctuations of the terms in expansion 
(25) for a hole in the t-J-model lead to poor statistics at large times [53], 
though, SO method is capable of recovering energy and Z-factor even from 
GF known only at small imaginary times [53]. 

2.3 Advantages and Drawbacks of DMC-SO Machinery 

Among numerical methods, capable of obtaining quantitative results in the 
problem of exciton (1) and (2), one can list time-dependent density func- 
tional theory [89], Hanke-Sham technique of correcting particle-hole excita- 
tion energy [90, 91], and approaches directly solving Bcthc-Salpetcr equation 
[92, 93, 94]. The latter ones provide rather accurate information on the two- 
particle GF. However, usage of finite mesh in direct /reciprocal space, which 
is avoided in DMC method, leads to its' failure in Wannier regime [93]. 

In contrast to DMC method, none of the traditional numeric methods 
can give reliable results for measurable properties of excited states of polaron 
at arbitrary range of electron-phonon interaction for the macroscopic system 
in the thermodynamic limit. Exact diagonalization method [95, 96, 97, 98] 
can study excited states though only on rather small finite size systems and 
results of this method are not even justified in the variational sense in the 
thermodynamic limit [99] . There is a batch of rather effective variational "ex- 
act translation" methods [99, 100, 101, 102, 103] where basis is chosen in 
the momentum space and, hence, the variational principle is applied in the 
thermodynamic limit. Although these methods can reveal few discrete excited 
states, its fail for long-range interaction and for dispersive, especially acoustic 
phonons due to catastrophic growth of variational basis. A non perturbative 
theory, which is able to give information about spectral properties in the ther- 
modynamic limit at least for one electron, is Dynamical Mean Field Theory 
[104, 105, 106, 107]. However it gives an exact solution only in the case of 
infinite dimension which does not correspond to a realistic system and can be 
considered only as a guide for extrapolation to finite dimensions [108]. 
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Recently developed cluster perturbation theory, where exact diagonaliza- 
tion of a cluster is further improved by taking into account inter-cluster inter- 
action [109, 110, 111, 112, 113], is applicable for study of the excited states, 
but limited to one-dimensional lattices or two-dimensional systems with short- 
range interaction. Traditional density-matrix renormalization group method 
[114, 115, 116, 117, 118] is very effective though mostly limited to one- 
dimensional systems and ladders. Finally, recently developed path integral 
quantum Monte Carlo algorithm [119, 120, 121, 122] is valid for any dimen- 
sion and properly takes into account quasi long-range interactions [123]. Path 
integral method is capable of obtaining the density of states [119, 120] and 
isotope exponents [121, 124]. However calculations of measurable characteris- 
tics of excited states, such as ARPES or optical conductivity, by this method 
were never reported. 

In conclusion, none of methods, except DMC-SO combination, can obtain 
at the moment approximation-free results for measurable physical quantities 
for a few QPs interacting with a macroscopic bosonic bath in the thermody- 
namic limit. Indeed, there are limitations of the DMC and SO methods. DMC 
method does not work in many-fermion systems due to sign problem and SO 
method fails at high temperatures, comparable to energies of dominant spec- 
tral peaks, because even very small statistical noise of GFs turns Fredholm 
equation (13) into essentially "ill defined" problem [84]. 

3 Spectral Properties of the Frohlich Polaron 

Before development of DMC-SO methods, the information on the excited 
states of polaron models, especially the Frohlich one, was very limited. Knowl- 
edge of LF was based on results of infinite-dimensions approximation [125], 
exact diagonalization [126, 96, 97, 97], or strong coupling expansion [127]. No 
one of the above techniques was capable of obtaining the LF of polaron with- 
out approximations, especially for long-range interaction where difficulties of 
traditional numerical methods dramatically increase. In a similar way, optical 
conductivity (OC) of Frohlich model was known only in strong coupling ex- 
pansion approximation [128], within the framework of the perturbation theory 
[129], or was based on the variational Feynman path integral technique [75]. 
In this sect, we consider exact DMC-SO results on LF [45] and OC [48, 57] of 
Frohlich polaron model. 

3.1 Lehmann Function of the Frohlich Polaron 

The perturbation theory expression for the high-energy part (oj > 0) of the 
LF for arbitrary interaction potential V(\ q |) reads [45] (frequency of the 
optical phonon w p /, is set to unity) 

i k=0 (^ > 0) = 1 2 VZ V T | V(V^T)) | 2 9(w 1) . (34) 
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Fig. 2. Comparison of the numerical results (solid lines) and the perturbation 
theory (dashed lines) for the LFs of the Frohlich model with a — 0.05 (a) and the 
short-range interaction model with a — 0.05 and n — 1 (b). LFs of Frohlich polaron 
for a = 0.5 (c), a = 1 (d) and a = 2 (e). Energy is measured from that of the 
ground state of the polaron. The initial fragment of the LF for a — 1 is shown in 
the inset (f). 



1/2 

i (2\/2a:7r) [q 2 + k 2 ) -1 / 2 , reducing to the Frohlich one when k — > 0, is 

L k=o (w<0) = ^^(5 f uj + a I ■ (35) 

Comparison of low-energy parts of the LF of the Frohlich model, obtained 
by DMC-SO and taken from (35), shows perfect agreement for a = 0.05: the 
accuracy for the polaron energy and Z-factor is about 10 -4 . On the other hand, 
high-energy part of numeric result (Fig. 2) significantly deviates from that of 
the analytic expression (35). This is not surprising since for Frohlich polaron 
the perturbation theory expression is diverging as uj — > cj p h and, therefore 
the perturbation theory breaks down. When perturbation theory is obviously 
valid, e.g. for the case of finite k = 1, there is a perfect agreement between 
analytic expression and DMC-SO results (Fig. 2b). Note that the high-energy 
part of Lk=o( ti; ) is successfully restored by SO method despite the fact that 
the total weight of the feature for a = 0.05 is less than 10~ 2 . 

The main deviation of the actual LF from the perturbation theory result 
is the extra broad peak in the actual LF at u> <~ 3.5. To study this feature 
Lk=o( w ) was calculated for a = 0.5, a = 1, and a — 2 (Fig. 2c-e). The peak 
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Fig. 3. Evolution of spectral density with a in the cross-over region from interme- 
diate to strong couplings. The polaron ground state peak is shown only for a = 8. 
Note that the spectral analysis still resolves it, despite its very small weight < 1CT 3 . 

is seen for higher values of the interaction constant and its weight grows with 
a. Near the threshold, w = 1, LF demonstrates the square-root dependence 
- (Fig. 2f). 

To trace the evolution of the peak at higher values of a the LF was cal- 
culated [45] for a = 4, a — 6, and a = 8 (Fig. 3). At a = 4 the peak at 
uj <~ 4 already dominates. Moreover, a distinct high-energy shoulder appears 
at a = 4, which transforms into a broad peak at u> ~ 8.5 in the LF for a = 6. 
The LF for a = 8 demonstrates further redistribution of the spectral weight 
between different maxima without significant shift of the peak positions. 

3.2 Optical Conductivity of the Frohlich Polaron: Validity of the 
Franck- Condon Principle in the Optical Spectroscopy 

The FC principle [130, 131] and its validity have been widely discussed in stud- 
ies of optical transitions in atoms, molecules [132, 133], and solids [134, 9]. 
Generally, the FC principle means that if only one of two coupled subsys- 
tems, e.g. an electronic subsystem, is affected by an external perturbation, 
the second subsystem, e.g., the lattice, is not fast enough to follow the re- 
construction of the electronic configuration. It is clear that the justification 
for the FC principle is the short characteristic time of the measurement pro- 
cess T mp <C T ic , where r mp is related to the energy gap between the initial 
and final states, AE, through the uncertainty principle: r mp ~ h/(AE) and 
Ti C is the time necessary to adjust the lattice when the electronic component 
is perturbed. Then, the spectroscopic response considerably depends on the 
value of the ratio T mp /Ti C For example, in mixed valence systems, where the 
ionic valence fluctuates between configurations / 5 and f 6 with characteristic 
time t^ ~ 10 _13 s, spectra of fast and slow experiments are dramatically dif- 
ferent [135, 136]. Photoemission experiments with short characteristic times 
T m p ~ 10~ 16 s (FC regime), reveal two lines, corresponding to f 5 and f 6 
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states. On the other hand slow Mossbauer isomer shift measurements with 
T mp ~ 10~ 9 s show a single broad peak with mean frequency lying between 
signals from pure f 5 and f 6 shells. Finally, according to paradigm of mea- 
surement process time, magnetic neutron scattering with r mp r ic revealed 
both coherent lines with all subsystems dynamically adjusted and broad inco- 
herent remnants of strongly damped excitation of / 5 and / 6 shells [137, 138]. 
Actually, the meaning of the times n c and r mp varies with the system and 
with the measurement process. 

To study the interplay between measurement process time T mp and ad- 
justment time T ic , the OC of the Frohlich polaron was studied in [57] from 
the weak to the strong coupling regime by three methods. DMC method gives 
numerically exact answer which is compared with memory function formalism 
(MFF), which is able to take dynamical lattice relaxation into account, and 
strong coupling expansion (SCE) which assumes FC approach. It was found 
that near critical coupling a c w 8.5 a dramatic change of the OC spectrum 
occurs: dominating peak of OC splits into two satellites. In this critical regime 
the upper (lower) one quickly decreases (increases) it's spectral weight as the 
value of coupling constant increases. Besides, while OC follows prediction of 
MFF at a < a c , its dependence switches to that predicted by SCE for larger 
couplings. It was concluded that, for the OC measurement of polaron, the 
adjustment time r ic w h/V is set by typical nonadiabatic energy V. Nona- 
diabaticity destroys FC classification at a < a c while FC principle rapidly 
regains its validity at large couplings due to fast growth of energy separation 
between initial and final states of optical transitions. 

Comparison of exact DMC-SO data for OC with existing results of ap- 
proximate methods showed [48] that the Feynman path integral technique 
[75] of Devreese, De Sitter, and Goovaerts, where OC is calculated starting 
from the Feynman variational model [139], is the only successfully describing 
evolution of the energy of the main peak in OC with coupling constant a (see 
[58]). However, starting from the intermediate coupling regime this approach 
fails to reproduce the peak width. Subsequently, the path integral approach 
was rewritten in terms of MFF [140]. Then, in [57] the extended MFF for- 
malism, which introduces dissipation processes fixed by exact sum rules, was 
developed [141]. 

As shown in Fig. 4a, in the weak coupling regime, the MFF, with or with- 
out dissipation, is in very good agreement with DMC data, showing significant 
improvement with respect to weak coupling perturbation approach [129] which 
provides a good description of OC spectra only for very small values of a. For 
1 < a < 8, where standard MFF fails to reproduce peak width (Fig. 4b-d) 
and even the peak position (Fig. 4c), the damping, introduced to extended 
MFF scheme, becomes crucial. Results of extended MFF are accurate for the 
peak energy and quite satisfactory for the peak width (Fig. 4b-e). Note that 
the broadening of the peak in DMC data is not a consequence of poor quality 
of analytic continuation procedure since DMC-SO methods is capable of re- 
vealing such fine features as 2- and 3-phonon thresholds of emission (Fig. 4b). 
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Fig. 4. Comparison of the optical conductivity calculated by DMC method (circles), 
extended MFF (solid line), and DSG [75, 140] (dashed line) for different values of 
a. The slanted arrows indicate 2- and 3-phonon thresholds of absorption. 
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Fig. 5. (a)-(c) Comparison of the optical conductivity calculated within the DMC 
method (circles), the extended MFF (solid line), and SCE (dotted line) for different 
values of a. (d) The energy of lower- and higher-frequency features (circles and tri- 
angles, respectively) compared with the FC transition energy with the SCE (dashed 
line) and with the energy of the peak obtained from the extended MFF (solid line) . 
In the inset, the weights of FC and adiabatically connected transitions are shown as 
a function of a (for r\ = 1.3.) 



However, a dramatic change of OC occurs around critical coupling strength 
a c « 8.5. The dominating peak of OC splits into two ones, the energy of lover 
one corresponding to the predictions of SCR expansion and that of upper 
one obeying extended MFF value (Fig. 5a). The shoulder, corresponding to 
dynamical extended MFF contribution, rapidly decreases it's intensity with 
increase of a and at large a (Fig. 5b-c) the OC is in good agreement with 
strong coupling expansion, assuming FC scheme. Finally, comparing energies 
of the peaks, obtained by DMC, extended MFF and FC strong coupling ex- 
pansion (Fig. 5d), we conclude that at critical coupling a c ~ 8.5 the spectral 
properties rapidly switch from dynamic, when lattice relaxes at transition, to 
FC regime, where nuclei are frozen in initial configuration. 

In order get an idea of the FC breakdown authors of [57] consider the fol- 
lowing arguments. The approximate adiabatic states are not exact eigenstates 
of the system. These states are mixed by nondiagonal matrix elements of the 
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nonadiabatic operator V and exact eigenstates are linear combinations of the 
adiabatic wavefunctions. Being interested in the properties of transition from 
ground (g) to an excited (ex) state, whose energy correspond to that of the OC 
peak, it is necessary to consider mixing of only these states and express exact 
wavefunctions as a linear combinations [142, 143] of ground and excited adi- 
abatic states. The coefficients of superposition are determined from standard 
techniques [142, 143] where nondiagonal matrix elements of the nonadiabatic 
operator [142] are expressed in terms of matrix elements of the kinetic energy 
operator M, the gap between excited and ground state AE = E ex — E g and 
the number rip of phonons in adiabatic state: 



The extent to which lattice can follow transition between electronic states, 
depends on the degree of mixing between initial and final exact eigenstates 
through the nonadiabatic interaction. If initial and ground states are strongly 
mixed, the adiabatic classification has no sense and, therefore, the FC pro- 
cesses have no place and lattice is adjusted to the change of electronic states 
during the transition. In the opposite limit adiabatic approximation is valid 
and FC processes dominate. The estimation for the weight of FC component 
Ifc [57] is equal to unity in the case of zero mixing and zero in the case 
of maximal mixing. The weight of adiabatically connected (AC) transition 
I Ac — 1 ~ Ifc is defined accordingly. Non-diagonal matrix element M is pro- 
portional to the root square of a with a coefficient 77 of the order of unity. 
In the strong coupling regime, assuming that AE w -fa 2 (7 k 0.1 from MC 
data), and rip w AE (rip 3> 1), one gets 



where T mp = l/AE and r ic — 1/D. For r\ of the order of unity one obtains 
qualitative description of a rather fast transition from AC- to FC-dominated 
transition, when Ipc and Iac exchange half of their weights in the range 
of a from 7 to 9. The physical reason for such quick change is the faster 
growth of energy separation AE <~ a 2 compared to that of the matrix ele- 
ment M ~ a 1 / 2 . Finally, for large couplings, initial and final states become 
adiabatically disconnected. The rapid AC-FC switch has nothing to do with 
the self-trapping phenomenon where crossing and hybridization of the ground 
and an excited states occurs. This phenomenon is a property of transition 
between different states and related to the choice whether lattice can or can 
not follow adiabatically the change of electronic state at the transition. 




(36) 



I FC = [1 + 4(t„ ip /t ic ) 2 ] 



(37) 



4 Self- Trapping 



In this section we consider the self-trapping (ST) phenomenon which, due to 
essential importance of many-particle interaction of QP with bosonic bath of 
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macroscopic system, was never addressed by exact method before. We start 
with a basic definition of the ST phenomenon and introduce the adopted 
criterion for it's existence. Then, generic features of ST are demonstrated 
on a simple model of Rashba-Pekar exciton-polaron in Sect. 4.1. It is shown 
in Sect. 4.2 that the criterion is not a dogma since even in one dimensional 
system, where ST is forbidden by criterion of existence, one can observe all 
main features of ST due to peculiar nature of electronic states. 

In general terms [7, 80], ST is a dramatic transformation of a QP properties 
when system parameters are slightly changed. The physical reason of ST is a 
quantum resonance, which happens at some critical interaction constant 7 C , 
between "trapped" (T) state of QP with strong lattice deformation around 
it and "free" (F) state. Naturally, ST transition is not abrupt because of 
nonadiabatic interaction between T and F states and all properties of the QP 
are analytic in 7 [144]. At small 7 < j c , ground state is an F state which is 
weakly coupled to phonons while excited states are T states and have a large 
lattice deformation. At critical couplings 7 -f c a crossover and hybridization 
of these states occurs. Then, for 7 > -f c the roles of the states exchange. The 
lowest state is a T state while the upper one is an F state. 

First, and up to now the only quantitative criterion for ST existence was 
given in terms of the ground state properties in the adiabatic approximation. 
This criterion considers stability of the delocalized state in undistorted lattice 
A = with respect to the energy gain due to lattice distortion A' 7^ 0. ST 
phenomenon occurs when completely delocalized state with A = is separated 
from distorted state with A' ^ by a barrier of adiabatic potential. One of 
these states is stable while another one is meta-stable. The criterion of barrier 
existence is defined in terms of the stability index 

s = d- 2(1 + 0, (38) 

where d is the system dimensionality. Index I determines the range of the force 
lim g ^o ipiq) ~ <l~ l , where tp(R) is the kernel of interaction U(R n ) = tp(R n — 
R n /)i/(R„/) connecting potential £/(R„) with generalized lattice distortion 
v(H n r ) [7]. The barrier exists for s > and does not exist for s < 0. The 
discontinuous change of the polaron state, i.e. ST, occurs in the former case 
while does not happen in the latter case. When s = 0, this scaling argument 
alone can not conclude the presence or absence of the ST and more detailed 
discussion for each model is needed. 

4.1 Typical Example of the Self- Trapping: Rasba-Pekar 
Exciton-Polaron 

Classical example of a system with ST phenomenon is the three dimensional 
continuous Rasba-Pekar exciton-polaron in the approximation of intraband 
scattering, i.e. when polar electron-phonon interaction (EPI) with dispersion- 
less optical phonons u> p h = 1 does not change the wave function of internal 




Fig. 6. The ground-state energy (a), effective mass (b), and average number of 
phonons as function of coupling constant (c). Partial weights of n-phonon states (d) 
in the polaron ground state (k = 0) at 7 = 18 (circles), 7 = 18.35 (squares), and 
7 = 19 (diamonds). Dotted line in panel (a) is the result of strong coupling limit 
and dashed line is the result of perturbation theory. 



electron-hole motion. System is defined as a structureless QP with dispersion 
e(k) = k 2 /2 and short range coupling to phonons [54, 7]. General criterion of 
the existence of ST is satisfied for three dimensional system with short range 
interaction [54, 7, 50] and, thus, one expects to observe typical features of the 
phenomenon. 

It is shown [54] that in the vicinity of the critical coupling j c ss 18 the 
average number of phonons (N) in (18) and effective mass m* quickly in- 
crease in the ground state by several orders of magnitude (Fig. 6b-c). Besides, 
a quantum resonance between polaronic phonon clouds of F and T state is 
demonstrated. Distribution of partial n-phonon contributions Z^ k=0 ' (n) in 
(17) has one maximum at n = in the weak coupling regime, which cor- 
responds to weak deformation, and one maximum at n > 1 in the strong 
coupling regime, which is the consequence of a strong lattice distortion. How- 
ever, due to F-T resonance there are two distinct peaks at n — and n » 1 
for 7 w 7 C (Fig. 6d). 

Near the critical coupling 7 C the LF of polaron has several stable states 
(Fig. 7 a-b) below the threshold of incoherent continuum E gs +cu p h- Any state 
above the threshold is unstable because emission of a phonon with transition 
to the ground state at k = with energy E gs is allowed. On the other hand, 
decay is forbidden by conservation laws for states below the threshold. De- 
pendence of the energies of ground and excited resonances on the interaction 
constant resembles a picture of crossing of several states interacting with each 
other (Fig. 7c). 

According to the general picture of the ST phenomenon, lowest F state in 
the weak coupling regime at k = has small effective mass m* w m of the 
order of the bare QP mass m. To the contrary, the effective mass of excited 
state m* > m is large. Hence, below the critical coupling the energy of the 
F state, which is lowest at k = 0, has to reach a flat band of T state at 
some momentum. Then, F and T state have to hybridize and exchange in 
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Fig. 7. LF Z/( k=0 )(cj) at critical coupling 7 = j c (a) and for 7 > 7c (b). Energy is 
counted from the polaron ground state, (c) Dependence of energy of ground state 
(squares) and stable excited states (circles, diamonds, and triangles) on the coupling 
constant. Dashed line is the threshold of the incoherent continuum. Dependence 
of energy (d) and average number of phonons (e) on the wave vector at 7 < j c 
(circles and rectangles). Dashed line is the effective mass approximation E ^ = 
E gs + k 2 /2m* for parameters E gs — —3.7946 and m* = 2.258, obtained by DMC 
estimators for given value of 7. Dotted line is a parabolic dispersion law which is 
fitted to last 4 points of energy dispersion curve with parameters E\ — —3.5273 and 
mj = 195. Empty square is the energy of first excited stable state at zero momentum 
obtained by SO method. 

energy. DMC data visualize this picture (Fig. 7 d-e). After F state crosses the 
flat band of excited T state, the average number of phonons increases and 
dispersion becomes flat. 

It is natural to assume that above the critical coupling the situation is 
opposite: ground state is the T state with large effective mass while excited 
F state has small, nearly bare, effective mass. Indeed, this assumption was 
confirmed in the framework of another model which is considered in Sect. 6.1. 
Moreover, it was shown that in the strong coupling regime excited resonance 
inherits not only bare effective mass around k = but the whole dispersion 
law of the bare QP [49]. 

4.2 Degeneracy Driven Self- Trapping 

According to the criterion (38), ST phenomenon in one-dimensional sys- 
tem does not occur. Although this statement is probably valid for the 
case of single band in relevant energy range, it is not the case for the 
generic multi-band cases. This fact has been unnoticed for many years, 
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which prevented the proper explanation of puzzling physics of quasi-one- 
dimensional compound Anthracenc-PMDA, although it's optical properties 
[65, 145, 146, 147, 66, 148] directly suggested resonance of T and F states. 
The reason is that in Anthracene-PMDA, in contrast to conditions at which 
criterion (38) is obtained, there are two, nearly degenerate exciton bands. 
Then, one can consider quasi-degenerate self-trapping mechanism when ST 
phenomenon is driven by nondiagonal interaction of phonons with quaside- 
generate exciton levels [52]. Such mechanism was already suggested for expla- 
nation of properties of mixed valence systems [143] though it's relevance was 
never proved by an exact approach. 




Fig. 8. Dependence of energy (a) and average number of phonons (b) on the non- 
diagonal coupling constant A12 at An = and A22 = 0.25. Phonon distributions in 
polaron cloud below ST point at A12 = 1.0125 (c), at ST point at A12 = 1.0435 (d), 
and above ST coupling at A12 = 1.0625 (e). 

The minimal model to demonstrate the mechanism of quasi-degenerate 
self-trapping involves one optical phonon branch with frequency ui p h = 0.1 
and two exciton branches with energies £1,2(9) = ^1,2 + 2[1 — cos(g)], where 
Ai = and A2 = 1. Presence of short range diagonal 722 and nondiagonal 
712 interactions (with corresponding dimensionless constants A22 = 722/(2^) 
and A12 = 7i2/(2w)) leads to classical self-trapping behavior even in one- 
dimensional system [52] (see Fig. 8). 

5 Exciton 

Despite numerous efforts over the years, there has been no rigorous tech- 
nique to solve for exciton properties even for the simplest model (l)-(2) which 
treats electron-electron interactions as a static renormalized Coulomb poten- 
tial with averaged dynamical screening. The only solvable cases are the Frenkel 
small-radius limit [67] and the Wannier large-radius limit [68] which describe 
molecular crystals and wide gap insulators with large dielectric constant, re- 
spectively. Meanwhile, even the accurate data for the limits of validity of the 
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Wannier and Frenkel approximations have not been available. As discussed in 
Sects. 1.2 and 2.3, semianalytic approaches has little to add to problem when 
quantitative results are needed whereas traditional numerical methods fail to 
reproduce them even in the Wannier regime. To the contrary, DMC results 
do not contain any approximation. 




Fig. 9. Panel (a): dependence of the exciton binding energy on the bandwidth 
E c = E v for conduction and valence bands. The dashed line corresponds to the 
Wannier model. The solid line is the cubic spline, the derivatives at the right and left 
ends being fixed by the Wannier limit and perturbation theory, respectively. Inset 
in panel (a): the initial part of the plot. Panel (b): the wave function of internal 
motion in real space for the optically forbidden monopolar exciton. Panels (c)-(e): 
the wave function of internal motion in real space: (c) Wannier [E c = E v = 60]; (d) 
intermediate [E c = E v — 10]; (e) near- Frenkel [E c — E v — 0.4] regimes. The solid 
line in the panel (c) is the Wannier model result while solid lines in other panels are 
to guide the eyes only. 

To study conditions of validity of limiting regimes by DMC method, 
electron-hole spectrum of three dimensional system was chosen in the form 
of symmetric valence and conduction bands with width E c and direct gap E g 
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at zero momentum [46]. For large ratio W — E c /E g , when W > 30, exci- 
ton binding energy is in good agreement with Wannier approximation results 
(Fig. 9a) and probability density of relative electron-hole motion corresponds 
(Fig. 9c) to hydrogen-like result. The striking result is the requirement of 
rather large valence and conduction bandwidths (W > 20) for applicability of 
Wannier approximation. For smaller values of W the binding energy and wave 
function of relative motion (Fig. 9d) deviate from large radius results. In the 
similar way, conditions of validity of Frenkcl approach are rather restricted 
too. Moreover, even strong localization of wave function does not guarantee 
good agreement between exact and Frenkel approximation result for binding 
energy. At 1 < W < 10 the wave function is already strongly localized though 
binding energy considerably differs from Frenkel approximation result. For 
example, at W = 0.4 relative motion is well localized (Fig. 9e) whereas the 
binding energy of Frenkel approximation is two times larger than exact result 
(Inset in Fig. 9a). 

A study of conditions necessary for formation of charge transfer exciton in 
three dimensional systems is crucial to finalize protracted discussion of numer- 
ous models concerning properties of mixed valence semiconductors [149]. A 
decade ago unusual properties of SmS and SmB 6 were explained by invoking 
the excitonic instability mechanism assuming charge-transfer nature of the 
optically forbidden exciton [150, 151]. Although this model explained quanti- 
tatively the phonon spectra [152, 153], optical properties [154, 155], and mag- 
netic neutron scattering data [138], it's basic assumption has been criticized 
as being groundless [156, 157]. To study excitonic wavefunction, dispersions 
of valence and conduction bands were chosen as it is typical for mixed valence 
materials: almost flat valence band is separated from broad conduction band, 
having maximum in the centre and minimum at the border of Brillouin zone 
[46]. Results presented in Fig. 9b support assumption of [150, 151] since wave 
function of relative motion has almost zero on-site component and maximal 
charge density at near neighbors. 

6 Polarons in Undoped High Temperature 
Superconductors 

It is now well established that the physics of high temperature superconduc- 
tors is that of hole doping a Mott insulator [158, 159, 160]. Even a single 
hole in a Mott insulator, i.e. a hole in an antiferromagnct in case of infinite 
Hubbard repulsion U, is substantially influenced by many-body effects [10] be- 
cause it's jump to a neighboring site disturbs antifcrromagnetic arrangement 
of spins. Hence a thorough understanding of the dynamics of doped holes in 
Mott insulators has attracted a great deal of recent interest. The two major 
interactions relevant to the electrons in solids are electron-electron interac- 
tions (EEI) and electron-phonon interactions (EPI). The importance of the 
former at low doping is no doubt essential since the Mott insulator is driven 
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by strong Hubbard repulsion, while the latter was considered to be largely 
irrelevant to superconductivity based on the observations of a small isotope 
effect on the optimal T c [161] and an absence of a phonon contribution to the 
resistivity (for review see [162]). 

On the other hand, there are now accumulating evidences that the EPI 
plays an important role in the physics of cuprates such as (i) an isotope effect 
on superfluid density p s and T c away from optimal doping [163], (ii) neutron 
and Raman scattering [164, 165, 166] experiments showing strong phonon soft- 
ening with both temperature and hole doping, indicating that EPI is strong 
[167, 168]. Furthermore, the recent studies of cuprates by the angle resolved 
photocmission spectroscopy (ARPES), which spectra are proportional to the 
LF (7) [32], resulted in the discovery of the dispersion "kinks" at around 40- 
70meV measured from the Fermi energy, in the correct range of the relevant 
oxygen related phonons [169, 170, 171]. These particular phonons - oxygen 
buckling and half-breathing modes are known to soften with doping [172, 164] 
and with temperature [170, 171, 172, 164, 165, 166] indicating strong cou- 
pling. The quick change of the velocity can be predicted by any interaction of 
a quasiparticle with a bosonic mode, either with a phonon [170, 171] or with 
a collective magnetic resonance mode [173, 174, 175]. However, the recently 
discovered "universality" of the kink energy for LSCO over the entire doping 
range [176] casts doubts on the validity of the latter scenario as the energy 
scale of the magnetic excitation changes strongly with doping. 

Besides, measured in undoped high T c materials ARPES revealed appar- 
ent contradiction between momentum dependence of the energy and linewidth 
of the QP peak. On the one hand the experimental energy dispersion of the 
broad peak in many underdoped compounds [31, 177] obeys the theoretical 
predictions [178, 179], whereas the experimental peak width is comparable 
with the bandwidth and orders of magnitude larger than that obtained from 
theory of Mott insulator [53]. Early attempts to interpret this anomalously 
short lifetime of a hole by an interaction with additional nonmagnetic bosonic 
excitations, e.g. phonons [180], faced generic question: is it possible that in- 
teraction with media leaves the energy dispersion absolutely unrenormalizcd, 
while, induces a decay which inverse life-time is comparable or even larger 
than the QP energy dispersion? A possibility of an extrinsic origin of this 
width can be ruled out since the doping induces further disorder, while a 
sharper peak is observed in the overdoped region. 

In order to understand whether phonons can be responsible for peculiar 
shape of the ARPES in the undoped cuprates, the LF of an interacting with 
phonons hole in Mott insulator was studied by DMC-SO [49]. The case of the 
LF of a single hole corresponds to the ARPES in an undoped compound. For 
a system with large Hubbard repulsion U, when U is much larger than the 
typical bandwidth W of noninteracting QP, the problem reduces to the t-J 
model [181, 182, 158, 11] 
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4-j = -t c l c js + ( s * s j - n ^ n j/ A ) ■ ( 39 ) 

(ij) 

Here Cj„ is projected (to avoid double occupancy) fermion annihilation op- 
erator, rii (< 2) is the occupation number, Si is spin 1/2 operator, J is an 
exchange integral, and (ij) denotes nearest-neighbor sites in two dimensional 
square lattice. Different theoretical approaches revealed [158, 183, 53] basic 
properties of the LF. The LF has a sharp peak in the low energy part of the 
spectrum which disperses with a bandwidth Wju ~ 2 J and, therefore, the 
large QP width in experiment can not be explained. More complicated tt't"-J 
model takes into account hoppings to the second t' and third t" nearest neigh- 
bors and, hence, dispersion of the hole changes [184, 185, 186, 178, 179, 32]. 
However, for parameters, which are necessary for description of dispersion 
in realistic high T c superconductors [31, 178], peak in the low energy part 
remains sharp and well defined for all momenta [187]. 

After expressing spin operators in terms of Holstcin-Primakoff spin wave 
operators and diagonalizing the spin part of Hamiltonian (39) by Fourier 
and Bogoliubov transformations [188, 10, 189, 190], tt't"-J Hamiltonian is 
reduced to the boson-holon model, where hole (annihilation operator is h^) 
with dispersion e(k) = At' cos(k x ) cos(fc y ) + 2t"[cos(2fc :E )+cos(2fc y )] propagates 
in the magnon (annihilation operator is au) bath 

Hi j = ]T e(k)ftt fc k + ]T wtatok (40) 



with magnon dispersion Wk = 2Jyl — 7^, where 7k = (cosk x + cosfc y )/2. 
The hole is scattered by magnons as described by 



tfh-m = i V-V2£ Mk;q 
k,q 



fr k frk-qak + h.c. (41) 



with the scattering vertex M k , q . Parameters t, t' and t" are hopping ampli- 
tudes to the first, second and third near neighbors, respectively. If hopping 
integrals t' and t" are set to zero and bare hole has no dispersion, the problem 
(40-41) corresponds to t-J model. 

Short range interaction of a hole with dispersionless optical phonons 
jje-ph — Ho^2 k fo^&k of the frequency S7o is introduced by Holstein Hamil- 
tonian 

' (42) 



k,q 



where a is the momentum and isotope independent coupling constant, M is 
the mass of the vibrating lattice ions, and J?o is the frequency of dispersionless 
phonon. The coefficient in front of square brackets is the standard Holstein in- 
teraction constant 7 = 0-/ ^J(2MQo). In the following we characterize strength 
of EPI in terms of dimcnsionless coupling constant A = 7 2 /4i/?o- Note, if in- 
teraction with magnetic subsystem (41) is neglected and hole dispersion e(k) 
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is chosen in the form e(k) = 2t[cos(k x ) +cos(k y )], the problem (40), (42) cor- 
responds to standard Holstein model where hole with near neighbor hopping 
amplitude t interacts with dispersionless phonons. 

We consider the evolution of ARPES of a single hole in t- J-Holstein model 
(40)-(42) from the weak to the strong coupling regime and dispersion of the LF 
in the strong coupling regime in Sect. 6.1. It occurs that properties of the LF in 
the strong coupling regime of the EPI explain the puzzle of broad lineshape 
in ARPES in underdoped high T c superconductors. Therefore, in order to 
suggest a crucial test for the mechanism of phonon-induced broadening, we 
present calculations of the effect of the isotope substitution on the ARPES in 
Sect. 6.2. 

6.1 Spectral Function of a Hole Interacting with Phonons in the 
t-J Model: Self- Trapping and Momentum Dependence 

Previously, the LF of t-J- Holstein model was studied by exact diagonalization 
method on small clusters [191] and in the non-crossing approximation (NCA) 6 
for both phonons and magnons [192, 193]. However, the small system size in 
exact diagonalization method implies a discrete spectrum and, therefore, the 
problem of lineshape could not be addressed. The latter method omits the 
FDs with mutual crossing of phonon propagators and, hence, is an invalid 
approximation for phonons in strong and intermediate couplings of EPI. This 
statement was demonstrated by DMC, which can sum all FDs for Holstein 
model both exactly and in the NCA [49] . Exact results and those of NCA are 
in good agreement for small values A < 0.4 and drastically different for A > 1. 
For example, for flo/t — 0.1 exact result shows a sharp crossover to strong 
coupling regime for A > AJ, k 1.2 whereas NCA result does not undergo such 
crossover even at A = 100. On the other hand, NCA is valid for interaction of a 
hole with magnons since spin S=l/2 can not flip more than once and number 
of magnons in the polaronic cloud can not be large. Note that the t-J- Holstein 
model is reduced to problem of polaron which interacts with several bosonic 
fields (3)- (4). 

DMC expansion in [49] takes into account mutual crossing of phonon prop- 
agators and, in the framework of partial NCA, neglects mutual crossing of 
magnon propagators, to avoid sign problem. NCA for magnons is justified for 
J/t < 0.4 by good agreement of results of NCA and exact diagonalization on 
small clusters [188, 10, 194, 195, 190]. Recently results of exact diagonalization 
were compared in the limit of small EPI for t- J-Holstein model, boson-holon 
model (40-42) without NCA, and boson-holon model with NCA [196]. Al- 
though agreement is not so good as for pure t-J model, it was concluded that 
NCA for magnons is still good enough to suggest that one can use NCA for a 
qualitative description of the t- J-Holstein model. 



6 NCA is equivalent to self-consistent Born approximation (SCBA) 
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Fig. 10. (a) The LF of a hole in the ground state k = (tt/2, n, 2) at J/t = 0.3 and 
A = 0. Low energy part of the LF of a hole in the ground state k = (n/2, tt, 2) at 
J/t = 0.3: (b) A = 0; (c) A = 0.3; (d) A = 0.4; (e) A = 0.46. Dependence on coupling 
strength A at J/t = 0.3: (f) energies of lowest LF resonances; (g) Z-factor of lowest 
peak; (h) average number of phonons (N). 

Figures lOa-e show low energy part of LF in the ground state at k = 
(n/2, tt/2) in the weak, intermediate, and strong coupling regimes of inter- 
action with phonons. Dependence on the coupling constant of energies of 
resonances (Fig. lOf), Z k= ( 7r / 2,7r ' /2 )-factor of lowest peak (Fig. lOg), and aver- 
age number of phonons in the polaronic cloud (N) (Fig. lOh) demonstrates a 
picture which is typical for ST (see [80, 54] and Sect. 4). Two states cross and 
hybridize in the vicinity of critical coupling constant AJ?_j s=s 0.38, Z k= ( 7r / 2 ' 7r / 2 )- 
factor of lowest resonance sharply drops and average number of phonons in 
polaronic cloud quickly rises. According to the general understanding of the 
ST phenomenon, above the critical couplings A > A£_j one expects that the 
lowest state is dispersionless while the upper one has small effective mass. 
This assumption is supported by the momentum dependence of the LF in 
the strong coupling regime (Fig. lla-e). Dispersion of upper broad shake-off 
Franck-Condon peak nearly perfectly obeys relation 

£k = e m in + W j/ t /5{[cos k x +cos ky} 2 + [cos(k x + ky)+cos(k x -ky)} 2 / 4}, (43) 

which describes dispersion of the pure t-J model in the broad range of ex- 
change constant 0.1 < J/t < 0.9 [194] (Fig. llf). Note that this property of 
the shake-off peak is general for the whole strong coupling regime (Fig. llf). 

Momentum dependence of the shake-off peak, reproducing that of the free 
particle, is the direct consequence of the adiabatic regime. Actually, phonon 
frequency J7 is much smaller than the coherent bandwidth 2J of the t-J 
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Fig. 11. The LF of a hole at J/t = 0.3 and A = 0.46: (a) full energy range 
for k = (7r/2, 7r/2); (b-e) low energy part for different momenta. Slanted arrows 
show broad peaks which can be interpreted in ARPES spectra as coherent (C) and 
incoherent (I) part. Vertical arrows in panels (b)-(e) indicate position of "invisible" 
lowest resonance, (f) Dispersion of resonances energies at J/t = 0.3: broad resonance 
(filled circles) and lowest polaron pole (filled squares) at A = 0.46; broad resonance 
(open circles) and lowest polaron pole (open squares) at A = 0.4. The solid curves 
are dispersions (43) of a hole in pure t-J model at J/t — 0.3 (Wj/t=o.a = 0.6): 
Emin = —2.396 (smin = —2.52) for dotted (solid) line. Panel (g) shows ground state 
potential Q 2 /2 (solid line), excited state potential without relaxation D + Q 2 /2 
(dashed line), and relaxed excited state potential D + (Q — A) 2 /2 — A 2 /2 (dotted 
line) . 



model, giving the adiabatic ratio Q/2J = 1/6 <C 1. Besides, as experience 
with the OC of the Frohlich polaron (Sect. 3.2) shows, there is one more 
important parameter in the strong coupling limit. Namely, the ratio between 
measurement process time T mp = Ti/AE where AE is the energy separation of 
shake-off hump from the ground state pole, and that of characteristic lattice 
time t w l/i?o is much less than unity Hence, fast photoemission probe sees 
the ions frozen in one of possible configurations [197]. The LF in the FC limit 
is a sum of transitions between a lower E\ ow (Q) and an upper E up (Q) sheets 
of adiabatic potential, weighted by the adiabatic wave function of the lower 
sheet | Viow(Q) | 2 [198]. If EPI is absent both in initial E low (Q) = Q 2 /2 
and final E up (Q) = V + Q 2 /2 states, the LF is peaked at the energy V. 
Then, if there is EPI AE up (Q) — ~\Q only in the final state, i.e. when hole 
is removed from the Mott insulator, the upper sheet of adiabatic potential 
E up (Q) = V- A 2 /2 + (Q- A) 2 /2 has the same energy V at Q = 0. Since the 
probability function | Viow(<3) | 2 has maximum at Q = 0, the peak of the LF 
broadens but it's energy does not shift [198] (Fig. llg). 
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Behavior of the LF is the same as observed in the ARPES of undoped 
cuprates. The LF consists of a broad peak and a high energy incoherent con- 
tinuum (see Fig. 11a). Besides, dispersion of the broad peak "c" in Figs. 11 
reproduces that of sharp peak in pure t-J model (Fig. llb-f). The lowest dis- 
persionless peak, corresponding to small radius polaron, has very small weight 
and, hence, can not be seen in experiment. On the other hand, according to ex- 
periment, momentum dependence of spectral weight Z^ of broad resonance 
exactly reproduces dispersion of Z^-factor of pure t-J model. The reason 
for such perfect mapping is that in adiabatic case £2 /2J <C 1 all weight of 
the sharp resonance in t-J model without EPI is transformed at strong EPI 
into the broad peak. This picture implies that the chemical potential in the 
heavily underdoped cuprates is not connected with the broad resonance but 
pinned to the real quasiparticle pole with small Z- factor. This conclusion was 
recently confirmed experimentally [177]. 

Comparing the critical EPI for a hole in the i-J-Holstcin model (40-42) 
AJ?_j i=a 0.38 and that for Holstein model w 1.2 with the same value of 
hopping t, we conclude that spin-hole interaction accelerates transition into 
the strong coupling regime. The reason for enhancement of the role of EPI is 
found in [196]. Comparison of the EPI driven renormalization of the effective 
mass in t-J- Holstein and Holstein model shows that large effective mass in the 
t-J model is responsible for this effect. The enhancement of the role of EPI 
by EEI takes place at least for a single hole at the bottom of the t-J band. 
Had the comparison been made with half-filled model, the result would have 
been smaller enhancement or no enhancement at all [199]. On the other hand, 
coupling constant of half-breathing phonon is increased by correlations [200] . 
Finally, we conclude that effect of enhancement of the effective EPI by EEI is 
not unambiguous and depends on details of interaction and filling. However, 
this effect is present for small filling in the t- J-Holstein model. 

6.2 Isotope Effect on ARPES in Underdoped High- Temperature 
Superconductors 

The magnetic resonance mode and the phonon modes are the two major 
candidates to explain the "kink" structure of the electron energy dispersion 
around 40-70 meV below the Fermi energy, and the isotope effect (IE) on 
ARPES should be the smoking-gun experiment to distinguish between these 
two. Gweon et al. [201] performed the ARPES experiment on 18 -replaced 
Bi2212 at optimal doping and found an appreciable IE, which however can 
not be explained within the conventional weak-coupling Migdal-Eliashberg 
theory. Namely the change of the spectral function due to 18 -replacement 
has been observed at higher energy region beyond the phonon energy (~ 
60meV). This is in sharp contrast to the weak coupling theory prediction, i.e., 
the IE should occur only near the phonon energy. Hence the IE in optimal 
Bi2212 remains still a puzzle. On the other hand, the ARPES in undoped 
materials, as described in Sect. 6.1, has recently been understood in terms of 
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the small polaron formation [49, 202, 198]. Therefore, it is essential to compare 
experiment in undoped systems with presented in this Sect. DMC-SO data, 
where theory can offer quantitative results. 

In addition to high-T c problem, strong EPI mechanism of ARPES spec- 
tra broadening was considered as one of alternative scenarios for diatomic 
molecules [203] , colossal magnetoresistive manganites [34] , quasi-one-dimensi- 
onal Peierls conductors [37, 38], and Verwey magnetites [39]. Therefore, exact 
analysis of the IE on ARPES at strong EPI is of general interest for conclusive 
experiments in a broad variety of compound classes. 

Dimcnsionlcss coupling constant A = 7 2 /4£i? in (42) is an invariant quan- 
tity for the simplest case of IE. Indeed, assuming natural relation Q ~ 1/\[M 
between phonon frequency and mass, we find that A does not depend on 
the isotope factor k; so = fi/flo = \/M a /M, which is defined as the ratio of 
phonon frequency in isotope substituted (Q) and normal (Oq) systems. We 
chose adopted parameters of the tt't"-J model which reproduce the experi- 
mental dispersion of ARPES [178]: J/t = 0.4, t'/t = -0.34, and t" /t = 0.23 . 
The frequency of the relevant phonon [32] is set to Q$/t = 0.2 and the isotope 
factor k; so = y/ 16/18 corresponds to substitution of O 18 isotope for O 16 . 

To sweep aside any doubts of possible instabilities of analytic continuation, 
we calculate the LF for normal compound (K n or = 1), isotope substituted 
(k; so = y/16/18) and "anti-isotope" substituted (/t an t = yl8/16) compounds. 
Monotonic dependence of LF on k ensures stability of analytic continuation 
and gives possibility to evaluate the error-bars of a quantity A using quantities 

*4iso -4nor, -^nor -^ant, and (.Also -^ant)/2. 

Since LF is sensitive to strengths of EPI only for low frequencies [55] , we 
concentrate on the low energy part of the spectrum. Figure 12 shows IE on the 
hole LF for different couplings in nodal and antinodal points, respectively. The 
general trend is a shift of all spectral features to larger energies with increase 
of the isotope mass (n < 1). One can also note that the shift of broad FCP is 
much larger than that of narrow real-QP peak. Moreover, for large couplings 
A the shift of QP energy approaches zero and only decrease of QP spectral 
weight Z is observed for larger isotope mass. On the other hand, the shift 
of FCP is not suppressed for larger couplings. Except for the LF in nodal 
point at A = 0.62 (Fig. 12a, b), where LF still has significant weight of QP 
6- functional peak, there is one more notable feature of the IE. With increase 
of the isotope mass the height of FCP increases. Taking into account the 
conservation law for LF L^(lu) — 1 and insensitivity of high energy part 
of LF to EPI strength [55] , the narrowing of the FCP for larger isotope mass 
can be concluded. To understand the trends of the IE in the strong coupling 
regime we analyze the exactly solvable independent oscillators model (IOM) 
[60]. The LF in IOM is the Poisson distribution 

L(oj) = expKoA] £ [ ^^GM - (44) 
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Fig. 12. Low energy part of hole LFs: normal compound (solid line), isotope sub- 
stituted compound (dotted line) and "antiisotope" substituted compound (dashed 
line). LFs at different couplings in the nodal (a, c, e) and antinodal (g, i, 1) points. 
Insets (b, d, f, h, k) show low energy peak of real QP. 

where £ = 7o/^o = 4tA/J? is dimensionless coupling constant for normal 
system and G K .i(<+>) = 5[cj + 4tA — S7 kI] is the S- function. The properties of the 
Poisson distribution quantitatively explain many features of the IE on LF 7 . 

The energy ujqp = — AtX of the zero-phonon line I = in (44) depends 
only on isotope independent quantities which explains very weak isotope de- 
pendence of QP peak energy in insets of Fig. 12. Besides, change of the zero- 
phonon line weight obeys relation Z^/Z^} T = exp [— £q(1 — K )/ K ] m 
IOM. These IOM estimates agree with DMC data within 15% in the nodal 
point and within 25% in the antinodal one. IE on FCP in the strong cou- 
pling regime follows from the properties of zero Mo = J^°° L(lo)(Lu = 1, first 

Mi = f*°° u)L(iv)cL = 0, and second M 2 = f^°° uj 2 L(uj)cLj = k^Hq mo- 
ments of shifted Poisson distribution (44). Moments Mq and Mi establish 
relation V = hf£ p /hlg p = \ « 1.03 between heights of FCP in normal 
and substituted compounds. DMC data in the antinodal point perfectly agree 
with the above estimate for all couplings. This is consistent with the idea that 
the anti-nodal region remains in the strong coupling regime even though the 
nodal region is in the crossover region. In the nodal point DMC data well 
agree with IOM estimate for A = 0.75 (V w 1.025) whereas at A = 0.69 and 

7 Cautions should be made about approximate form of EPI (42) . Strictly speaking, 
actual momentum dependence of the interaction constant a [204, 205] can slightly 
change the obtained differences between nodal and antinodal points though the 
general trends have to be left intact because ST is caused solely by the short 
range part of EPI [80]. 
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Fig. 13. (a) Energies of ground state and broad peaks for normal (triangles), 
isotope substituted (circles) and "antiisotopc" substituted (diamonds) compounds. 
Comparison of IOM estimates (lines) with DMC data in the nodal (squares) and 
antinodal (diamonds) points: (b) shift of the FCP top, (c) FCP leading edge at 1/2 
of height, and (d) FCP leading edge at 1/3 of height. 



A = 0.62 influence of the ST point leads to anomalous values of V: V 1.07 
and V w 0.98, respectively. Shift of the low energy edge at half maximum 
A\/2 must be proportional to change of the root square of second moment 
Z\ y rpj = Vlo^o[l — \fR\. As we found in numeric simulations of (44) with 
Gaussian functions 8 Q k .i{lo), relation A X j 2 ~ A^j^/2 is accurate to 10% for 
0.62 < A < 0.75. Also, simulations show that the shift of the edge at one 
third of maximum /A1/3 obeys relation Ai/ 3 w Ayjjj- DMC data with IOM 
estimates are in good agreement for strong EPI A = 0.75 (Fig. 13). How- 
ever, shift of the FCP top A p and Ai/ 2 are considerably enhanced in the 
self-trapping (ST) transition region. The physical reason for enhancement of 
IE in this region is a general property regardless of the QP dispersion, range 
of EPI, etc. The influence of nonadiabatic matrix element, mixing excited 
and ground states, on the energies of resonances essentially depends on the 
phonon frequency. While in the adiabatic approximation ST transition is sud- 
den and nonanalytic in A [80] , nonadiabatic matrix elements turn it to smooth 
crossover [144]. Thus, as illustrated in Fig. 13a, the smaller the frequency the 
sharper the kink in the dependence of excited state energy on the interaction 
constant 

In the undoped case the present results can be directly compared with 
the experiments. It is found that the IE on the ARPES lincshapc of a sin- 
gle hole is anomalously enhanced in the intermediate coupling regime while 
can be described by the simple independent oscillators model in the strong 
coupling regime. The shift of FCP top and change of the FCP height are rele- 
vant quantities to pursue experimentally in the intermediate coupling regime 
since IE on these characteristics is enhanced near the self trapping point. In 



Results are almost independent on the parameter r\ of the Gaussian distribution 
Q K ,i(u>) = l/(r]V2n) exp(-[w + 4f A - Q kI]/(2t] 2 )) in the range [0.12, 0.2]. 



Spectroscopic Properties of Polarons by Exact Monte Carlo 35 

contrast, shift of the leading edge of the broad peak is the relevant quantity 
in the strong coupling regime since this value increases with coupling as VA- 
These conclusions, depending on the fact whether self trapping phenomenon 
is encountered in specific case, can be applied fully or partially to another 
compounds with strong EPI [34, 37, 38, 39]. 

6.3 Conclusions and Perspectives 

In this article, we have focused mainly on the polaron problem in strongly 
correlated systems. This offers an approach from the limit of low carrier con- 
centration doped into the (Mott) insulator, which is complementary to the 
conventional Eliashbcrg-Migdal approach for the EPI in metals. In the latter 
case, we have the Fermi energy sp as a relevant energy scale, which is usually 
much larger than the phonon frequency ft®. In this case, the adiabatic Migdal 
approximation is valid and the vertex corrections, which correspond to the 
multi-phonon cloud and are essential to the self-trapping phenomenon, are 
suppressed by the ratio Qq/ef- Therefore an important issue is the crossover 
from the strong coupling polaronic picture to the weak coupling Eliashberg- 
Migdal picture. This occurs as one increases the carrier doping into the insu- 
lator. As is observed by ARPES experiments in high temperature supercon- 
ductors, the polaronic states continue to survive even at finite doping [177]. 
This suggests a novel polaronic metallic state in underdoped cuprates, which 
is common also in CMR manganites [36] and is most probably universal in 
transition metal oxides. In the optimal and overdoped region, the Eliashberg- 
Migdal picture becomes appropriate [170, 171], but still a nontrivial feature of 
the EPI is its strong momentum dependence leading to the dichotomy between 
the nodal and anti-nodal regions. It is an interesting observation that the high- 
est superconducting transition temperature is attained at the crossover region 
between the two pictures above, which suggests that both the itinerancy and 
strong coupling to the phonons are essential to the quantum coherence. It 
should be noted that this crossover occurs in a nontrivial way also in the mo- 
mentum space, i.e., the nodal and anti- nodal regions behave quite differently 
as discussed in Sect. 6.2. However, the relevance of the EPI to the high T c 
superconductivity is still left for future investigations. 

We hope that this article convinces the readers the vital role of ARPES 
experiments and numerically exact solutions to the EPI problem, the com- 
bination of them offers a powerful tool for the momentum-energy resolved 
analysis of these rather complicated strongly correlated electronic systems. 
This will pave a new path to the deeper understanding of the many-body 
electronic systems. 
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